// Copyright 2007 by Rory Plaire (codekaizen@gmail.com)

using System;
using System.Collections.Generic;
using System.Text;

namespace Novas.Net
{
	public static class Computations
	{
		public static double[] Precession(double firstEpoch, double[] position, double secondEpoch)
		{
			double xx, yx, zx, xy, yy, zy, xz, yz, zz, t, t1, t02, t2, t3,
			   zeta0, zee, theta, cz0, sz0, ct, st, cz, sz;

			/*
			   't' and 't1' below correspond to Lieske's "big T" and "little t".
			*/

			t = (firstEpoch - Constants.Epoch) / 36525.0;
			t1 = (secondEpoch - firstEpoch) / 36525.0;
			t02 = t * t;
			t2 = t1 * t1;
			t3 = t2 * t1;

			/*
			   'zeta0', 'zee', 'theta' below correspond to Lieske's "zeta-sub-a",
			   "z-sub-a", and "theta-sub-a".
			*/

			zeta0 = (2306.2181 + 1.39656 * t - 0.000139 * t02) * t1
				  + (0.30188 - 0.000344 * t) * t2 + 0.017998 * t3;

			zee = (2306.2181 + 1.39656 * t - 0.000139 * t02) * t1
				+ (1.09468 + 0.000066 * t) * t2 + 0.018203 * t3;

			theta = (2004.3109 - 0.85330 * t - 0.000217 * t02) * t1
				  + (-0.42665 - 0.000217 * t) * t2 - 0.041833 * t3;

			zeta0 /= Constants.RadiansToSeconds;
			zee /= Constants.RadiansToSeconds;
			theta /= Constants.RadiansToSeconds;

			/*
			   Precalculate trig terms.
			*/

			cz0 = Math.Cos(zeta0);
			sz0 = Math.Sin(zeta0);
			ct = Math.Cos(theta);
			st = Math.Sin(theta);
			cz = Math.Cos(zee);
			sz = Math.Sin(zee);

			/*
			   Precession rotation matrix follows.
			*/

			xx = cz0 * ct * cz - sz0 * sz;
			yx = -sz0 * ct * cz - cz0 * sz;
			zx = -st * cz;
			xy = cz0 * ct * sz + sz0 * cz;
			yy = -sz0 * ct * sz + cz0 * cz;
			zy = -st * sz;
			xz = cz0 * st;
			yz = -sz0 * st;
			zz = ct;

			/*
			   Perform rotation.
			*/

			double[] result = new double[3];

			result[0] = xx * position[0] + yx * position[1] + zx * position[2];
			result[1] = xy * position[0] + yy * position[1] + zy * position[2];
			result[2] = xz * position[0] + yz * position[1] + zz * position[2];

			return result;
		}
	}
}
